library(tidybulk)
library(tidySummarizedExperiment)
library(ggpubr)
library(GEOquery)
library(readxl)
library(tidyverse)
source('00_util_scripts/mod_bulk.R')
source('00_util_scripts/mod_bplot.R')

dna.rep <- read_xlsx('mission/HMCES_HIV/DNA Damage Repair related genes.xlsx',
                     col_names = 'gene') |>
  pull(gene)

# GC DZ enriched genes -----------
## GSE139891 basso20 ------------
read_delim('mission/HMCES_HIV/basso20/GSM4560816_GC3a_Ab_counts.txt.gz')
library(data.table)
library(Seurat)
source('00_util_scripts/mod_seurat.R')

ba20.path <- list.files('mission/HMCES_HIV/basso20/', full.names = T)

ba20.mtx <- ba20.path |>
  str_extract('(GC|DZ|LZ)\\d') |>
  set_names(x = ba20.path, nm = _) |>
  map(fread, .progress = T)

gene3 <- ba20.mtx[[3]]$V1 |>
  str_remove('.+;')

ens3 <- ba20.mtx[[3]]$V1 |>
  str_extract('ENSG\\d+')


tidy.basso <- function(mtx, name){
  mtx |>
    mutate(V1 = str_remove(V1, '.+;')) |>
    distinct(V1, .keep_all = T) |>
    column_to_rownames('V1') |>
    add_name_field(name)
}

ba20.rnm <- ba20.mtx |> 
  imap(tidy.basso, .progress = T)

ba20.rnm[[1]][1:4,1:4]['WASH7P',]

shared.gene <- ba20.rnm |>
  map(rownames) |>
  purrr::reduce(intersect)

ba20.shared <- ba20.rnm |>
  map(\(x)x[shared.gene,])

ba20.dgcm <- ba20.shared |>
  map(\(x)x |> as.matrix() |> as('dgCMatrix'))

sobj <- ba20.dgcm[[1]] |>
  RowMergeSparseMatrices(ba20.dgcm[-1]) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj %<>% filter(mito.ratio < 10)

sobj %<>% mutate(tissue = str_remove(orig.ident, '\\d'))

sobj %<>% NormalizeData()

sobj |> write_rds('mission/HMCES_HIV/basso20/basso20.rds')

dz.marker <- sobj |>
  FindMarkers(group.by = 'tissue',
              ident.1 = 'DZ') |>
  as_tibble(rownames = 'gene')

dz.marker |>
  filter(p_val_adj < .05, avg_log2FC > 0) |>
  write_csv('mission/HMCES_HIV/basso20/gc.dz.marker.csv')

dz.marker <-
  read_csv('mission/HMCES_HIV/basso20/gc.dz.marker.csv')

# bradley2018 --------
bra_meta <- read_delim('mission/HMCES_HIV/bradley2018.meta.txt',
                       col_names = c('acc','sample')) |>
  separate(sample, into = c('group','name'))

bradley <- read_csv('mission/HMCES_HIV/bradley2018_GSE115449_EffCountFile.csv.gz')

bradley <- bradley |>
  pivot_longer(2:last_col())

bradley <- bra_meta |>
  select(-acc) |>
  right_join(bradley)

bradley <- bradley$Transcript |>
  clusterProfiler::bitr(fromType = 'REFSEQ',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  dplyr::rename(Transcript = REFSEQ) |>
  left_join(bradley) |>
  as_tibble()

tdb_bra <- bradley |>
  tidybulk(.sample = name, .transcript = SYMBOL,
           .abundance = value)

## or use ncbi aligned counts
brad.nccount <- download_ncbi_counts('GSE115449')

tdb_bra <- brad.nccount |>
  pivot_longer(-1, names_to = 'acc') |>
  left_join(bra_meta) |>
  select(-acc) |>
  tidybulk(.sample = name, .transcript = Symbol,
           .abundance = value)

tdb_bra <- tdb_bra |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_bra |>
  plot_qc_bulk(value_scaled, group)
  
test_bra <- tdb_bra |>
  test_differential_abundance(~ 0 + group,
                              contrasts = c('groupBNab-groupControl'),
                              omit_contrast_in_colnames = T)

bra_sig <- test_bra |>
  keep_abundant() |>
  pivot_transcript()

bra_sig <- bra_sig |>
  mutate(SYMBOL = Symbol)

bra_sig |>
  write_csv('mission/HMCES_HIV/bradley.hiv.bnab.deg.csv')

bra_sig <-
  read_csv('mission/HMCES_HIV/bradley.hiv.bnab.deg.csv')

bra_sig <- bra_sig |>
  filter(SYMBOL %in% dz.marker$gene)

bra.dna <- bra_sig |>
  filter(SYMBOL %in% dna.rep) |>
  mutate(SYMBOL = fct_reorder(SYMBOL, logFC),
         type = case_when(SYMBOL == 'HMCES' ~ 'HMCES',
                          PValue < .05 & logFC > 0 ~ 'Upregulated',
                          PValue > .05 ~ 'NS',
                          .default = 'Downregulated'))

bra.dna |>
  filter(type != 'NS') |>
  arrange(desc(logFC)) |>
  mutate(rank = seq_along(SYMBOL)) |>
  filter(SYMBOL == 'HMCES') |>
  select(rank)

## barplot ------------
bra.dna |>
  filter(type != 'NS') |>
  ggplot(aes(logFC, SYMBOL, fill = type)) +
  geom_col() +
  theme_pubr(legend = 'right') +
  labs(title = 'DNA damage repair related genes in B cells of HIV patients\nbNAb vs control',
       x = 'log2FC', y = 'Gene', subtitle = 'GSE115449') +
  scale_fill_manual(values = c('royalblue','cyan4','red2')) +
  theme(axis.text.y = element_blank())

publish_pdf('mission/HMCES_HIV/figures/hiv.dna.repair.barplot.pdf',
            width = 140,100)

bra.dna |>
  filter(type == 'Upregulated') |>
  slice_max(logFC, n = 20) |>
  ggplot(aes(logFC, SYMBOL, fill = type)) +
  geom_col() +
  theme_pubr() +
  labs(title = 'Top 10 DNA mutation & repair genes in B cells of HIV patients\nbNAb vs control',
       x = 'log2FC', y = 'Gene', subtitle = 'GSE115449') +
  scale_fill_manual(values = c('royalblue','grey','red2'))

bra_sig |>
  filter(SYMBOL %in% dna.rep) |>
  mutate(SYMBOL = fct_reorder(SYMBOL, logFC),
         type = case_when(PValue < .05 & logFC > 0 ~ 'Upregulated',
                          PValue > .05 ~ 'NS',
                          .default = 'Downregulated')) |>
  filter(type != 'NS') |>
  ggplot(aes(logFC, SYMBOL, fill = type)) +
  geom_col() +
  theme_pubr() +
  labs(title = 'DNA mutation & repair genes in B cells of HIV patients\n bNAb vs control',
       x = 'log2FC', y = 'Gene', subtitle = 'GSE115449') +
  scale_fill_manual(values = c('royalblue','grey','red2'))

bra_sig |>
  filter(SYMBOL %in% dna.rep) |>
  select(SYMBOL, logFC, PValue) |>
  write_csv('mission/HMCES_HIV/hiv.GSE115449.dna.repair.barplot.csv')

bra.rep <-
  read_csv('mission/HMCES_HIV/hiv.GSE115449.dna.repair.barplot.csv')

bra.rep |>
  left_join(dz.marker, by = join_by(SYMBOL == gene)) |>
  filter(pct.1 > .5, avg_log2FC > .1, p_val_adj < .05,
         logFC > .1, PValue < .05)

tdb_bra |>
  filter(Symbol %in% dna.rep) |>
  select(Symbol, name, group, value_scaled) |>
  write_csv('hiv.GSE115449.dna.repair.rna_expr.csv')

bra_sig |>
  filter(SYMBOL == 'HMCES')

## volcano --------------
seek_name <- bra_sig |>
  mutate(gene = SYMBOL, avg_log2FC = logFC, p_val_adj = PValue) |>
  filter(SYMBOL %in% dz.marker$gene, SYMBOL %in% dna.rep, PValue < .05)

fc_thres <- .1
pval_thres <- .05

bra_sig |>
  filter(SYMBOL %in% dz.marker$gene) |>
  mutate(gene = SYMBOL, avg_log2FC = logFC, p_val_adj = PValue) |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj))) +
  geom_point(size = .5, alpha = .3, color = 'grey') +
  geom_point(data = seek_name, color = 'orange') +
  theme_pubr() +
  theme(legend.position = 'top', plot.margin = margin(), legend.box.margin = margin()) +
  #expand_limits(x = c(-symm_x_lim,symm_x_lim)) +
  labs(title = 'B cells of HIV patients with bNAbs vs control',
       subtitle = 'GSE115449')

publish_pdf('mission/HMCES_HIV/figures/hiv.dna.repair.volcano.png',
            width = 200, height = 150)

publish_pdf('mission/HMCES_HIV/figures/hiv.dna.repair.volcano.pdf',
            width = 200, height = 150)

tdb_bra |>
  filter(SYMBOL == 'HMCES') |>
  mutate(group = fct_relevel(group, 'Control')) |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3000, yend = 3000) +
  annotate(geom = 'text',x = 1.5, y=3150, size =5,
           label = 'logFC=0.135, P.adj=0.0734') +
  scale_color_discrete(label = c('Ctrl (n=46)', 'BNab (n=46)')) +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Bradley 2018 GSE115449')

tdb_bra |>
  filter(SYMBOL %in% dna.rep) |>
  mutate(group = fct_relevel(group, 'Control')) |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3000, yend = 3000) +
  annotate(geom = 'text',x = 1.5, y=3150, size =5,
           label = 'logFC=0.135, P.adj=0.0734') +
  scale_color_discrete(label = c('Ctrl (n=46)', 'BNab (n=46)')) +
  theme_pubr() +
  facet_wrap(~SYMBOL, scales = 'free_y') +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Bradley 2018 GSE115449')

## export scaled counts -------
candid <- c('HMCES','YWHAB','TDP2')

tdb_bra |>
  filter(SYMBOL %in% candid) |>
  select(SYMBOL, group, name, value_scaled) |>
  pivot_wider(names_from = SYMBOL, values_from = value_scaled) |>
  write_csv('mission/HMCES_HIV/hiv.GSE115449.3gene.expr.csv')

tdb_bra |>
  filter(SYMBOL %in% candid) |>
  mutate(logcount = log1p(value_scaled)) |>
  ggplot(aes(group, logcount, color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(vars(SYMBOL), scales = 'free_y')

# austin2019 -------
austin <- read_delim('mission/HMCES_HIV/Austin2019_GSE119234_LN_BcellSubset_raw_counts.txt.gz')

tdb_aus <- austin |>
  pivot_longer(2:last_col()) |>
  rename(symbol = ...1) |>
  mutate(sample = name) |>
  separate_wider_delim(name, delim = ':',
                       names = c('individual','cell_type')) |>
  tidybulk(.sample = sample,
           .transcript = symbol,
           .abundance = value)

tdb_aus <- tdb_aus |>
  identify_abundant(factor_of_interest = individual) |>
  scale_abundance()

tdb_aus |>
  ggplot(aes(value_scaled, group = sample, color = individual)) +
  geom_density() +
  scale_x_log10()

tdb_aus <- tdb_aus |>
  mutate(cell_type = str_remove(cell_type, '_HIV'))

tdb_aus <- tdb_aus |>
  mutate(infect = str_extract(individual, 'HD|HIV'))

tdb_aus |>
  reduce_dimensions(method = 'PCA') |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2, color = cell_type, shape = infect)) +
  geom_point()

tdb_aus |>
  filter(symbol == 'C3orf37') |>
  ggplot(aes(individual, cell_type, fill = log1p(value_scaled))) +
  geom_raster() +
  scale_fill_gradient(high = 'red', low = 'blue') +
  theme_pubr(legend = 'right') +
  labs(fill = 'log-Normalized expression',
       title = 'HMCES expression in B cell subsets',
       subtitle = 'Austin 2019 GSE119234')

tdb_aus |>
  filter(symbol == 'C3orf37') |>
  as_tibble() |>
  tidyHeatmap::heatmap(individual, cell_type, value_scaled,transform = log1p,
                       scale = 'column',
                       palette_value = c('blue','grey','red'))

## compare HIV GC vs Ctrl GC --------
aus_hivhd <- tdb_aus |>
  filter(cell_type == 'GC') |>
  test_differential_abundance(~0+infect,
                              contrasts = 'infectHIV-infectHD')

aus_hivhd |>
  pivot_transcript() |>
  filter(symbol == 'C3orf37') |>
  select(`FDR___infectHIV-infectHD`)

aus_hivhd |>
  filter(symbol == 'C3orf37') |>
  ggplot(aes(infect, value_scaled, color = infect)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .01) +
  annotate(geom = 'segment',x = 1, xend = 2, y=6100, yend =6100) +
  annotate(geom = 'text',x = 1.5, y=6300, size =5,
           label = 'logFC=-0.213, P.adj=0.65') +
  theme_pubr() +
  labs(y = 'Normalized expression', color = 'group',
       title = 'HMCES expression in GC B cell',
       subtitle = 'Austin 2019 GSE119234') +
  scale_color_manual(values = c('blue','red'),
                     label = c('HD (n=5)','HIV (n=4)'))

aus_hivhd |>
  filter(symbol == 'C3orf37' & infect == 'HIV')

## BNab-hi vs BNab-lo -----------
aus_bnab <- tdb_aus |>
  mutate(bnab = case_when(individual %in% c('HIV1','HIV4') ~ 'high',
                          infect == 'HIV' ~ 'low',
                          .default = NA)) |>
  filter(cell_type == 'GC' & infect == 'HIV') |>
  test_differential_abundance(~0+bnab,
                              contrasts = 'bnabhigh-bnablow')

aus_bnab |>
  pivot_transcript() |>
  filter(symbol == 'C3orf37') |>
  select(`logFC___bnabhigh-bnablow`,`FDR___bnabhigh-bnablow`)

aus_bnab |>
  filter(symbol == 'C3orf37') |>
  mutate(bnab = fct_relevel(bnab, 'low')) |>
  ggplot(aes(bnab, value_scaled, color = bnab)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .01) +
  annotate(geom = 'segment',x = 1, xend = 2, y=6150, yend =6150) +
  annotate(geom = 'text',x = 1.5, y=6350, size =5,
           label = 'logFC=0.505, P.adj=0.992') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in GC B cell',
       subtitle = 'Austin 2019 GSE119234') +
  scale_color_manual(values = c('red','blue'))

# martin2020 ---------
meta_mart <- getGEO('GSE141498')

meta_mart <-
meta_mart$GSE141498_series_matrix.txt.gz@phenoData@data |>
  select(title, `cohort:ch1`) |>
  as_tibble()

meta_mart <- meta_mart |>
  mutate(name = str_replace(title, ' for Patient ', '_'),
         id = str_extract(name, 'EMG.+'),
         cell = str_extract(name, '.+(?=_EMG)'),
         cohort = `cohort:ch1`,.keep = 'none')

martin <- read_csv('mission/HMCES_HIV/martin2020_GSE141498_RawCounts.csv.gz')

m <- martin |>
  separate_wider_delim(ID, names = c('ENSEMBL','SYMBOL'),delim = '|') |>
  select(-ENSEMBL) |>
  pivot_longer(-1) |>
  left_join(meta_mart)

tdb_mar <- m |>
  tidybulk(.sample = name,
           .transcript = SYMBOL,
           .abundance = value)

tdb_mar <- tdb_mar |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = cohort) |>
  scale_abundance()

tdb_mar |>
  ggplot(aes(value_scaled, group = id, color = cohort)) +
  geom_density() +
  scale_x_log10()

tdb_mar |>
  reduce_dimensions(method = 'PCA') |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2, color = cohort, shape = cell)) +
  geom_point()

tdb_mar |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(cell, id, fill = log1p(value_scaled))) +
  geom_raster() +
  scale_fill_gradient(high = 'red', low = 'blue') +
  theme_pubr(legend = 'right') +
  labs(fill = 'Normalized expression')

## compare Bnab vs no_Bnab --------
mar_neu <- tdb_mar |>
  filter(cell == 'Bcells') |>
  test_differential_abundance(~0+cohort,
                              contrasts = 'cohortNeut-cohortNonNeut')

mar_neu |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  select(`logFC___cohortNeut-cohortNonNeut`,`FDR___cohortNeut-cohortNonNeut`)


tdb_mar |>
  filter(cell == 'monocytes') |>
  test_differential_abundance(~0+cohort,
                              contrasts = 'cohortNeut-cohortNonNeut') |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  select(`logFC___cohortNeut-cohortNonNeut`,`FDR___cohortNeut-cohortNonNeut`)

mar_neu |>
  filter(SYMBOL == 'HMCES') |>
  mutate(cohort = fct_relevel(cohort, 'NonNeut')) |>
  ggplot(aes(cohort, value_scaled, color = cohort)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  annotate(geom = 'text',x = 1.5, y=800, size =5,
           label = 'logFC=-0.498, P.adj=0.721') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in B cell',
       subtitle = 'Martin 2020 GSE141498') +
  scale_color_manual(values = c('red','blue'),
                     label = c('No-bnAb (n=35)','bnAb (n=14)'))

tdb_mar |>
  filter(SYMBOL == 'HMCES' & cell != 'Bcells') |>
  mutate(cohort = fct_relevel(cohort, 'NonNeut')) |>
  ggplot(aes(cohort, value_scaled, color = cohort)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  annotate(geom = 'text',x = 1.5, y=800, size =5,
           label = 'NS') +
  theme_pubr() +
  facet_wrap(~cell) +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in other immune cells',
       subtitle = 'Martin 2020 GSE141498') +
  scale_color_manual(values = c('red','blue'),
                     label = c('No-bnAb (n=35)','bnAb (n=14)'))

# luo2022 ----------
luo22 <- read_excel('mission/HMCES_HIV/luo2022_GSE157198_readcount_GEO.xlsx')

luo_meta <- read_delim('mission/HMCES_HIV/luo2022_meta.txt',
                       col_names = c('acc','sample'))

luo_meta <- luo_meta |>
  mutate(name = str_extract(sample, 'S\\d+'),
         group = str_extract(sample, 'EC|Healthy|ART-na|ART') |>
           make.names(), .keep = 'none')

tdb_luo <- luo22$geneID |>
  clusterProfiler::bitr(fromType = 'ENSEMBL',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  rename(geneID = ENSEMBL) |>
  left_join(luo22) |>
  select(-geneID) |>
  as_tibble() |>
  pivot_longer(-1) |>
  left_join(luo_meta)

# identify abundance transcripts is highly recommended before scaling
tdb_luo <- tdb_luo |>
  tidybulk(.sample = name,
           .transcript = SYMBOL,
           .abundance = value) |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_luo |>
  reduce_dimensions(method = 'pca') |>
  pivot_sample() |>
  ggplot(aes(PC1, PC2, color = group)) +
  geom_point()

tdb_luo |>
  mutate(group = make.names(group)) |>
  test_differential_abundance(~0+group,
                              contrasts = 'groupEC-groupART.na') |>
  filter(SYMBOL == 'HMCES') |>
  as.data.frame()

tdb_luo |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  #annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  #annotate(geom = 'text',x = 1.5, y=800, size =5,
  #         label = 'NS') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in monocytes',
       subtitle = 'Luo 2022 GSE157198')

# parker2023 ----------
parker <- data.table::fread('mission/HMCES_HIV/parker2023_GSE199911_processeddata_genecounts.csv.gz')

parker[1:5,1:8]

parker <- parker |>
  mutate(group = ifelse(studygroup == 1, 'HIV', 'HC')) |>
  select(-c(V1,studygroup:age)) |>
  pivot_longer(where(is.numeric)) |>
  mutate(ENSEMBL = str_remove(name, '\\..+')) |>
  select(-name)

parker <- parker |>
  ensembl_to_symbol(.ensembl = ENSEMBL) |>
  filter(!is.na(transcript))

tdb_par <- parker |>
  as_tibble() |>
  tidybulk(.sample = sample,
           .transcript = SYMBOL,
           .abundance = value) |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_par |>
  reduce_dimensions(method = 'pca') |>
  pivot_sample() |>
  ggplot(aes(PC1, PC2, color = group)) +
  geom_point()

tdb_par |>
  test_differential_abundance(~0+group,
                              contrasts = 'groupHIV-groupHC') |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  as.data.frame()

tdb_par |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(group, value_scaled, color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3100, yend = 3100) +
  annotate(geom = 'text',x = 1.5, y=3200, size =5,
           label = 'NS') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Parker 2023 GSE199911')

# Wu-Chen 2021 COVID19 ---------
wu.salmon <- list.files('~/append-ssd/nextflowing/', recursive = T, full.names = T,
           pattern = 'salmon.merged.gene_counts_length') |>
  str_subset('batch.+tsv') |>
  map(read_delim)

wu.salmon <- wu.salmon |>
  purrr::reduce(left_join)

wu.meta <- read_delim('Archive/covid19/data/wu-chen-samples.txt')

wu.tidy <- wu.salmon |>
  select(-gene_id) |>
  pivot_longer(-1) |>
  mutate(id1 = str_extract(name, 'P\\d+')) |>
  left_join(wu.meta[c('id1','clinic')]) |>
  filter(!is.na(clinic))

wu.tidy |>
  write_csv('Archive/covid19/data/wu-chen-salmon.count.csv')

wu.tidy <-
  read_csv('Archive/covid19/data/wu-chen-salmon.count.csv')

wu.tdb <- wu.tidy |>
  select(-name) |>
  tidybulk(.sample = id1, .transcript = gene_name, .abundance = value)

wu.tdb <- wu.tdb |>
  preproc_bulk(clinic)

wu.tdb |>
  plot_qc_bulk(scaled_abundance = value_scaled, group = clinic)

## seems to have 2 batches independent of stages
wu.pca <- wu.tdb |>
  reduce_dimensions(method="PCA", .dims = 2) |>
  pivot_sample()

wu.kmeans <- wu.pca |>
  pull(PC2) |>
  kmeans(centers = 2)

wu.kmeans$cluster |>
  table()

wu.pca |>
  mutate(kmeans = wu.kmeans$cluster) |>
  ggplot(aes(PC1, PC2, color = clinic, shape = as.character(kmeans))) +
  geom_point()

wu.tdb <- wu.pca |>
  mutate(kmeans.pc2 = wu.kmeans$cluster) |>
  select(id1, kmeans.pc2) |>
  right_join(wu.tdb)

## test asympotomatic vs critical group
wu.res <- wu.tdb |>
  test_differential_abundance(~ kmeans.pc2 + clinic,
                              contrasts = 'clinic4_critical',
                              .sample = id1, .transcript = gene_name,
                              .abundance = value, omit_contrast_in_colnames = T)

wu.deg <- wu.res |>
  pivot_transcript(.transcript = gene_name) 

wu.deg |>
  mutate(gene = gene_name, avg_log2FC = -logFC, p_val_adj = FDR, .keep = 'none') |>
  filter(!is.na(avg_log2FC)) |>
  write_csv('mission/HMCES_HIV/covid19.ERP127339.DNA.repair.volcano.csv')

wu.deg <-
  read_csv('mission/HMCES_HIV/covid19.ERP127339.DNA.repair.volcano.csv')

## volcano ------------
seek_name <- wu.deg |>
  filter(gene %in% dz.marker$gene, gene %in% dna.rep, p_val_adj < .05)

fc_thres <- 1
pval_thres <- .05

wu.deg |>
  filter(gene %in% dz.marker$gene) |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj))) +
  geom_point(size = .5, alpha = .3, color = 'grey') +
  geom_point(data = seek_name, color = 'orange') +
  theme_pubr() +
  theme(legend.position = 'top', plot.margin = margin(), legend.box.margin = margin()) +
  #expand_limits(x = c(-symm_x_lim,symm_x_lim)) +
  labs(title = 'PBMC of COVID-19 patients of severity groups',
       subtitle = 'ERP127339')

publish_pdf('mission/HMCES_HIV/figures/covid19.dna.repair.volcano.png',
            width = 200, height = 150)

publish_pdf('mission/HMCES_HIV/figures/covid19.dna.repair.volcano.pdf',
            width = 200, height = 150)

## rank barplot -----------
wu.rep <- wu.deg |>
  filter(gene %in% dna.rep) |>
  mutate(gene = fct_reorder(gene, avg_log2FC),
         type = case_when(gene == 'HMCES' ~ 'HMCES',
                          p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj > .05 ~ 'NS',
                          .default = 'Downregulated')) |>
  write_csv('mission/HMCES_HIV/covid19.ERP127339.DNA.repair.barplot.csv')

wu.rep |>
  filter(type != 'NS') |>
  ggplot(aes(avg_log2FC, gene, fill = type)) +
  geom_col() +
  theme_pubr(legend = 'right') +
  labs(title = 'DNA damage repair genes expression in PBMC:\nAsymptomatic vs Critical',
       x = 'log2FC', y = 'Gene', subtitle = 'ERP127339') +
  scale_fill_manual(values = c('royalblue','cyan4','red2')) +
  theme(axis.text.y = element_blank())

publish_pdf('mission/HMCES_HIV/figures/covid19.dna.repair.barplot.pdf',
            width = 140,100)

wu.rep |>
  filter(type != 'NS') |>
  slice_max(avg_log2FC, n = 10) |>
  ggplot(aes(avg_log2FC, gene, fill = type)) +
  geom_col() +
  theme_pubr() +
  labs(title = 'Top 10 DNA damage repair genes expression\nPBMC of Asymptomatic vs Critical COVID-19',
       x = 'log2FC', y = 'Gene', subtitle = 'ERP127339') +
  scale_fill_manual(values = c('cyan4','red2'))

publish_pdf('mission/HMCES_HIV/figures/covid19.dna.repair.barplot.top10.pdf',
            width = 140,100)

## export 3 genes ----------
wu.cpm <- read_delim('Archive/covid19/results/wu-chen-cpm.tsv')

wu.cpm |>
  filter(SYMBOL %in% candid) |>
  pivot_longer(-1, names_to = 'run_accession') |>
  left_join(wu.meta) |>
  select(SYMBOL, id1, clinic, value) |>
  pivot_wider(names_from = SYMBOL, values_from = value) |>
  write_csv('mission/HMCES_HIV/covid19.ERP127339.3gene.expr.csv')

## heatmap of GOI ------
wu.asym.up <- wu.deg |>
  filter(avg_log2FC > 0, p_val_adj < .05)

wu.mean.expr <- wu.tdb |>
  filter(gene_name %in% dna.rep, gene_name %in% dz.marker$gene,
         gene_name %in% wu.asym.up$gene) |>
  summarise(mean_expr = mean(value_scaled), .by = c(gene_name, clinic)) |>
  as_tibble() 

wu.htmp <- wu.mean.expr |>
  tidyHeatmap::heatmap(gene_name, clinic, mean_expr, scale = 'row',
                       cluster_columns = F,
                       column_title = 'COVID-19 patients of severity groups',
                       row_title = 'Gene', name = 'z-score',
                       palette_value = c('blue','white','red'))

wu.htmp |>
  tidyHeatmap::save_pdf('mission/HMCES_HIV/figures/covid19.intersect.up10.heatmap.pdf',
                        width = 100, height = 100, units = 'mm')

wu.mean.expr |>
  mutate(zscore = scale(mean_expr)[,1], .by = gene_name) |>
  select(-mean_expr) |>
  pivot_wider(values_from = zscore, names_from = clinic) |>
  write_csv('covid19.ERP127339.DNA.repair.heatmap.csv')

## Venn of 4 plots --------
library(ggvenn)

dz.marker <-
  read_csv('mission/HMCES_HIV/basso20/gc.dz.marker.csv') |>
  filter()

dz.marker |>
  filter(gene == 'TDP2')

bra_sig <-
  read_csv('mission/HMCES_HIV/bradley.hiv.bnab.deg.csv') |>
  filter(logFC > 0, PValue < .05)

wu.deg <-
  read_csv('mission/HMCES_HIV/covid19.ERP127339.DNA.repair.volcano.csv') |>
  filter(avg_log2FC > 0, p_val_adj < .05)

venn.gene <- list('DZ GC marker' = dz.marker$gene,
     'COVID-19 asympotomatic' = wu.deg$gene,
     'DNA damage repair' = dna.rep,
     'HIV BNAb' = bra_sig$SYMBOL)

venn.gene |>
  ggvenn(set_name_size = 5, digits = 2)

publish_pdf('mission/HMCES_HIV/figures/venn4.pdf',
            width = 150, height = 150)

candid2 <- venn.gene |>
  purrr::reduce(intersect)

bra_sig |>
  filter(SYMBOL %in% candid2)

dz.marker |>
  filter(gene %in% dna.rep) |>
  write_csv('GC-DZ-DNA-repair-gene.csv')
